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Abstract Using data from nine microsatellite DNA loci and a population genetic approach, we evaluate the barrier 
effect of the Qinghai-Tibet Railway on toad-headed lizard, Phrynocephalus vlangalii. The study area is along a 20 
km stretch of the railway on northern Qinghai-Tibet Plateau, and this section of the railway was constructed between 
1958-1979. Both assignment tests and analysis of molecular variance (AMOVA) were used for data analysis. We found 
significant genetic differentiation between the populations from the study area and those from a further southeastern 
area, which are separated by a 20 km gap. This suggests the existence of population substructure at a fine-scale. 
However, we did not detect any difference between samples from the western and eastern sides of the railway within 
the study area, and concluded that the railway may not impose a significant barrier effect on these lizard populations 
at the present time. Available suitable habitat alongside the railway and bridge underpasses may have facilitated the 
gene exchange between the sides. The relatively short time since the completion of the railway may not allow the 
differentiation to accumulate to a detectable level. Since the Qinghai-Tibet Plateau maintains a unique and fragile 


ecosystem, long-term monitoring of such man-made landscape features is imperative for protecting this ecosystem. 
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1. Introduction consequences (Epps et al., 2005; Riley et al., 2006). 


Such impacts are particularly acute in a sensitive 
ecosystem, such as the Qinghai-Tibet Plateau. The recent 
construction of the Qinghai-Tibet Railway has created a 
storm of controversy surrounding its potential negative 


Roads and railways are among the most conspicuous 
man-made landscape features and they impose major 
ecological impacts on the ecosystem (Forman and 
Alexander, 1998; Trombulak and Frissell, 2000). Various 
effects have been proposed and examined, including 
large numbers of fatalities (Gibbs and Shriver, 2005), 
altering individual movement patterns and abundance 
(Shine et al., 2004; Fahrig and Rytwinski, 2009), and 
more importantly the “barrier effect” that subdivides 
populations with significant demographic and genetic 


impacts on the local ecosystem (Yin et al., 2006; Xia 
et al., 2007). Measures, such as wildlife passages, have 
been put into place to reduce such impacts and some of 
them appeared to be functioning (Yang and Xia, 2008). 
The long-term effects, however, are largely unknown. 
Although the completion of the railway was in 2006, the 
first phase construction of the railway (Xining to Golmud 
section) started in 1958 and was completed in 1979 
(Figure 1). This section of the railway has been in use for 
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more than 30 years and provides a unique opportunity to 
examine its impact on the local ecosystem in a relatively 
large time span. 

The Qinghai toad-headed lizard, Phrynocephalus 
vlangalii, is one of the dominant terrestrial species in 
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the Qinghai-Tibet Plateau. Its preferred habitat is sand 
dunes with loose substrate and little vegetation. Both 
males and females were observed to excavate burrows. 
Using eleven microsatellite DNA loci, Wang et al. (2009) 
reported a surprisingly large amount of genetic diversity 
and high population differentiation within this species. 
Furthermore, the species demonstrated strong territoriality 
(Qi et al., 2011). Such properties provide the potential 
to detect variations at a small spatial and temporal scale. 
Therefore, the species may serve as an excellent model 
system to examine the barrier effect of the Qinghai-Tibet 
Railway. 

Recent developments in individual assignment 
methods have also greatly improved the ability to detect 
population structure at a fine scale (Manel et al., 2003; 
Crosby et al., 2009). Using individual, rather than 
population, as the unit of examination, these methods 
do not require a pre-defined population boundary, which 
is particularly suitable for terrestrial species with low 
population density. Furthermore, these methods are very 
sensitive and are capable of detecting first generation 
migrants (Pritchard et al., 2000; Piry et al., 2004). 


Figure 1 The Xining-Golmud section of the Qinghai-Tibet 
Railway. Arrow indicates sand dunes, which are suitable habitat for 
toad-headed lizards. 


In this study, we analyzed data from nine microsatellite 
DNA loci using assignment tests to examine the barrier 
effect of the Qinghai-Tibet Railway on the toad-headed 
lizard, P. vlangalii. Specifically, we addressed the 
question whether the railway formed a genetic barrier to 
the gene flow of this lizard. If the Railway significantly 
impeded gene flow, we would expect the populations 
from both sides of the railway to have a different genetic 
make-up. On the other hand, if the railway did not have 


significant impact, the individuals from both sides would 
make up one well-mixed population. 


2. Materials and Methods 


2.1 Sampling A total of 170 individuals of P. vlangalii 
were collected along a 20 km stretch of the Qinghai- 
Tibet Railway in July 2010, including 91 from the 
western side and 79 from the eastern side of the railway. 
The sampling area is located in the Qaidam Basin with 
elevations between 2700-2800 m above sea level. The 
individuals were collected from a total of 34 sites, with 
1-11 individuals from each site. Our sampling strategy 
was to sample many sites but collect a few individuals 
from majority sites, which best suits the low population 
density nature of this lizard and the individual assignment 
analytic methods. 

We collected an additional 56 individual lizards from 
14 sites at an area located further southeast of our railway 
sampling area. There is a 20-25 km distribution gap 
between the railway sites and the southeastern sites; in the 
gap, the habitat is not suitable and lizards are rare. This 
additional collection is to confirm the natural occurrence 
of detectable population structure at a fine scale and to 
provide a comparison. All sampling sites and related 
information are presented in Figure 2 and Appendix I. 

All individuals were collected and euthanized in 
the field. Liver tissues were collected and preserved in 
95% ethanol. All voucher specimens are stored in the 


Figure 2 Sampling sites. Map is generated from Google maps 
using a satellite image. Circles show the sampling areas. Sites 1-35 
are nearby the railway and from both sides; actual sites are too close 
together to be accurately marked on this map. Sites 36—47 are for 
the purpose of comparison. Some sites are off the map. 
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2.2 Laboratory protocols Genomic DNA was extracted 
using TIANamp Genomic DNA kits (Tiangen Biotech) 
following the manufacturer’s recommendation. Nine 
polymorphic microsatellite DNA loci (Phr27, Phr75, 
Phr79, Phr63, Phr160, Pvms12, Pvms18, Pvms35, 
Pvms38) were amplified by polymerase chain reaction 
(PCR) using primers developed by Urquhart et al. (2005) 
and Zhan and Fu (2009). The forward primers were 
labeled with one of FAM, TAMRA or HEX florescence 
dyes. 

All PCR was conducted with individual pairs of 
primers in a 25 uL mix, including 20-50 ng of template 
DNA, 1U Taq DNA polymerase (TaKaRa), 2 mM MgCl, 
0.2 mM of each dNTP, and 10 pmol of each primer in 
1 x PCR buffer (TaKaRa). PCR cycling parameters were 
5 min for initial denaturation at 95°C, followed by 30 
cycles at 94°C denaturating for 30 sec, at 58°C annealing 
for 30 sec, at 72°C extension for 30 sec, and a final 
additional extension step at 72°C for 5 min. PCR products 
with different fluorescent labels were then pooled and 
genotyped on an ABI 3730 genetic analyzer (Applied 
Biosystems). Electropherograms were visualized and 
analyzed by Gene Marker version 1.6 (SoftGenetics). 

After the initial completion of genotyping all 
individuals, we randomly selected 20% of the samples 
and repeated the genotyping of these samples as a quality 
control measure. For three of the nine loci examined, 
a greater than 8% of inconsistency rate was detected, 
and therefore, all samples for these three loci were re- 
genotyped and inconsistent results were resolved. For the 
other six loci, the inconsistency rates were relatively low 
and the inconsistent results within these 20% repeated 
samples were resolved, but other samples were not re- 
examined. 


2.3 Data analysis Genetic diversity was assessed with 
three indices. The number of alleles (A), the observed 
heterozygosity (Ho) and the expected heterozygosity 
(H,,) were calculated using Genepop web version 4.0 
(Raymond and Rousset, 1995; Rousset, 2008). Departure 
from Hardy-Weinberg equilibrium (HWE) and linkage 
equilibrium (LE) was also examined using Genepop. 

To assess the potential barrier effect of the railway, two 
sets of assignment tests were conducted. Using Structure 
version 2.3.1 (Pritchard et al., 2009), a Bayesian model- 
based assignment test was employed to detect individual 
clustering and assign individuals to groups. The test uses 
multi-locus genotype data and a pre-defined number 


of clusters (K) to generate groupings that minimize the 
deviation from HWE and LE (Pritchard et al., 2000). 
Individuals are assigned probabilistically to one or 
more clusters based on their genotypic data. Two sets 
of Structure analyses were conducted: 1) all samples 
were included; and 2) only the samples from the railway 
sites (sites 1-34) were included. The admixture model 
and independent allele frequencies model were used, 
and a total of 50 independent runs with 300 000 burn- 
in iterations and 100 000 post burn-in iterations were 
conducted at each value of K. The ten highest LnP(D) 
values from the 50 runs were averaged to estimate 
the best K. Since we were only interested in detecting 
potential differentiation between the two sides of the 
railway, we restricted the range of K from one to four. 
We assessed both LnP(D) values and the individual 
assignment probabilities to obtain the best value of K. In 
general, higher LnP(D) and large increase in the LnP(D) 
values between successive changes in K indicate a better 
fit to the data. Furthermore, the best K value should 
produce an assignment plot within which the majority 
of individuals are assigned to a proper cluster with high 
probabilities. Admixed individuals with less than 80% 
probabilities should be minimal. 

A second set of assignment tests were conducted 
using GeneClass2 version 2.0.g (Piry et al., 2004). 
This method estimates the likelihood that a genotype 
originates from a particular population given the allele 
frequencies of the population, and is considered more 
sensitive than the method implemented in Structure 
(Waser and Hadfield, 2011). Only the samples from 
the railway sites (sites 1-34) were subjected to this 
analysis and individuals were pre-defined as from the 
“western population” or “eastern population” based on 
where they were collected. Rannala and Mountain’s 
(1997) resampling algorithm and Paetkau et al.’s (2004) 
simulation method (10 000 replicates) were used. 

A locus by locus analysis of molecular variance 
(AMOVA) was also conducted to assess the barrier effect 
of the railway using Arleoquin version 3.5 (Excoffier and 
Lischer, 2011). Permutation tests (10 000 replicates) were 
carried out at four hierarchical levels: among groups, 
within groups among sites, among individuals within 
sites, and within individuals. Similar to the Structure 
analysis, two sets of analyses were conducted. First, all 
samples were included to test the natural occurrence of 
population genetic structure at a fine scale. Two groups 
were defined, one including the railway sites (1-34) 
and the other including the southeastern sites (35—47). 
Second, only the samples from the railway sites (1-34) 
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were included to test whether the railway represents a 
genetic barrier. Again, two groups were defined: samples 
from the western side of the railway were defined as one 
and those from the eastern side were defined as the other. 


3. Results 


The number of alleles (4), the observed heterozygosity 
(Ho) and the expected heterozygosity (H,,) are presented 
in Table 1. In concordance with previous studies, we 
found relatively high genetic diversity within this species. 
For example, the number of alleles on each locus ranges 
from 18 to 33. 

Only the samples from the railway sites were subjected 
to HWE and LE tests. Individuals were pooled as a 
single population. Six loci showed significant departure 
from HWE after Bonferroni correction (Table 1). Two 
(out of 36) pairs of loci were in significant linkage 
disequilibrium. These results were expected because of 
Wahlund effect. Such departure provides information 
for the assignment tests. Previous studies of the same 
loci for P. vlangalii did not detect any consistent linage 
disequilibrium between populations (Wang et al., 2009), 
and therefore, we assumed these loci are not in physical 
linkage. 

The Structure assignment tests revealed two clusters 
when all samples were included in the analysis; however, 
one cluster was detected when only the samples from 
the railway sites were included. Figure 3 presents a 
plot of LnP(D) against K values of 1 to 4, and Figure 4 
presents the individual assignment probability bar graphs. 
When all samples were included, the LnP(D) showed 
a major increase when K changed from 1 to 2, and the 
increases of LnP(D) were minimal with the subsequent 
changes of K (Figure 3 A). In the bar graph, with K = 2, 
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Figure 3 Average likelihood values for K = 1—4 from Structure 
(ver. 2.3.1). LnP(D) values are average of the 10 highest likelihood 
values of 50 independent runs. A: All samples are included; B: Only 
the samples from the railway sites are included. 


approximately half of the individuals received assignment 
probabilities greater than 0.90, and three quarters of the 
individuals had probabilities greater than 0.80 (Figure 
4 A). Both measurements showed that K = 2 is the best 
fit for the data. It was also clear that the majority of 
individuals from the railway site made up one cluster 
and the majority from the southeastern sites formed the 
second cluster (Figure 4 A). When only the samples from 
the railway sites were included, the changes in K values 
did not produce a corresponding increase in LnP(D) 
values (Figure 3 B), suggesting that K = | is the best fit of 
data. In the bar graph, approximately 20% of individuals 
received assignment probabilities greater than 0.80, and 
no one received a probability of greater than 0.90 (Figure 
4 B). Considering both the change in LnP(D) values and 
the individual assignment probabilities, there is clear 
genetic differentiation between the railway sites and the 
southeastern sites, but all samples from the railway sites 
appeared to belong to one population. 

The results from the GeneClass2 analysis were 
consistent with these from the Structure analysis. Every 
individual from the railway sites received similar 


Table 1 Measurements of genetic diversity at nine microsatellite DNA loci. 


Populations Phr27 Phr75 Phr79 Phr63 Pyms12 Pyms18 Pyms35 Pvms38 Phr160 
Railway A 21 35 22 16 20 22. 19 31 16 
(west+east) Ho 0.659 0.800 0.765 0.694 0.612 0.553 0.618 0.894 0.806 
Hg 0.907 0.915 0.918 0.877 0.701 0.888 0.855 0.893 0.827 
Puy 0.0000 0.0000 0.0002 0.0000 0.0565 0.0000 0.0000 0.5834 0.2032 
Southeastern A 19 26 18 16 19 21 18 25 13 
Ho 0.491 0.527 0.673 0.782 0.655 0.582 0.800 0.927 0.691 
Ay 0.930 0.902 0.903 0.883 0.900 0.921 0.924 0.947 0.813 
A by locus 22 37 23 19 23 26 24 33 18 


A: Number of alleles; Ho: Observed heterozygosity; H,: Expected heterozygosity; P,,: P value of Hardy-Weinberg equilibrium test. Bold 
indicates significance after sequential Bonferroni correction at 0.01 level. 
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Figure 4 Individual assignment probability bar plots from Structure (ver. 2.3.1) for K = 2. Each vertical bar represents one individual. 
Individuals in both western and eastern sides of the railway are arranged in the order from north to south, and individuals from the 
southeastern sites are arranged in the order from west to east. Western = samples collected from the western side of the railway; Eastern 
= samples collected from the eastern side of the railway; Southeastern = samples collected from the southeastern sites. A: All samples are 


included; B: Only the samples from the railway sites are included. 


probabilities of being resident of both western and eastern 
populations; not a single individual received a particularly 
high probability of belonging to one side (probability 
difference > 0.5). This suggests that all samples from the 
railway sites make up one population. Twelve individuals 
were unlikely residents of either side the railway; their 
probabilities of being residents of any side were less than 
0.01. These individuals are likely migrants or descendants 
of recent migrants from other areas. 

AMOVA revealed a similar pattern (Table 2). When 
all samples were included, there was a small (3.27%) but 
significant difference between the railway site samples 
and the southeastern samples (P < 0.001; Table 2). When 
only the samples from the railway sites were included, 
there was no significant difference between the eastern 
sites and western sites (P > 0.874). Variation at all other 
hierarchical levels was significant. 


4. Discussion 


4.1 High genetic diversity This species has very high 
genetic diversity. Within the 20 km stretch of sampling 
area and 170 samples, the numbers of alleles per locus 
ranged from 16 to 35, with an average of 22.4 alleles 
per locus (Table 1). The diversity is much higher than 


that of many other lizard species in similar studies using 
microsatellite DNA loci. For example, the grand skink 
(Oligosoma grande) has an average number of alleles 
of 15.3 per locus (Berry et al., 2003). Ernst et al. (2004) 
found 8-20 alleles per locus in the Florida scrub lizard 
(Sceloporus woodi), and Urquhart et al. (2005) found 


Table 2 Results of the analysis of molecular variance (AMOVA). 


ete Sum of Variance Percentage 
Source of variation vee 
square components _ variation 
A: WHE vs. SE 
Among groups 28.600 0.13342 3.27337 
Among sites within groups 249.172 0.09868 2.42114 
Among individuals within sites 803.283 0.67201 16.48717" 
Within individuals 716.500 3.17183 77.81832" 
Total 1797.555 4.07594 
B: W vs. E 
Among groups 4.752 -0.00716 -0.18397 
Among sites within groups 177.624 0.10550 2.70947" 
Among individuals within sites 592.744 0.59535 15.29010° 
Within individuals 544.000 3.20000 82.18440" 
Total 1319.121 3.89368 


Grouping schemes are as following: A: All samples with two 
groups: 1) both sides of the railway sites, and 2) southeastern sites; 
B: Only the samples from the railway sites with two groups: 1) sites 
from western side of the railway, and 2) sites from eastern side of 
the railway. Average over 9 loci; number of permutations = 10 000; 
allowed level of missing data = 0.05. *: P < 0.001. 
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7-25 alleles per locus for a different toad-headed lizard 
species, P. przewalskii. The revealed high diversity in this 
species is consistent with a previous study (Wang et al., 
2009). They recorded an average of 32.5 alleles per locus 
over a much large sampling area. 

There is clear genetic structure among populations 
even within relatively short distance. The short gap (20— 
25 km) between the railway sites and the southeastern 
sites produced a detectable differentiation, by both the 
assignment tests and AMOVA. This observation also 
suggests that these tests are sensitive and able to pick up 
the signal from these genetic differences. 


4.2 Does the railway function as genetic barrier to 
the lizard populations? Our data suggest that at the 
present time the Qinghai-Tibet Railway does not impose 
a detectable “barrier effect” on the toad-headed lizards. 
Individuals from both sides of the railway share similar 
genetic make-up and form one population. There is 
no detectable sub-structure or clustering within the 
population. There are several possible explanations. 

The Qinghai-Tibet Railway may not significantly 
reduce the crossing of these lizards despite roads often 
obstructing reptile movement (Shine et al., 2004; 
Andrews and Gibbons, 2005). Railroads in particular 
have also been found to restrict gene flow in many other 
species (Bartoszek and Greenwald, 2009). Previous 
studies on large mammals found that the Qinghai-Tibet 
Railway indeed impeded crossing of these animals. For 
example, Yin et al. (2006) found that the Tibetan antelope 
(Pantholops hodgsoni), Tibetan gazelle (Procapra 
picticaudata) and Kiang (Equus kiang) were sensitive 
to roads, as well as other man-made structures. The 
raising embankment of the railway and human introduced 
objects, such as construction materials and the railway 
tracks, appeared to be the main cause of animal avoidance 
for large mammals (Xia et al., 2007). Such “alien” objects 
might not deter lizards as much as mammals. In addition, 
we noticed several bridges and culverts along the stretch 
of our sampling area, and lizards might use them as 
dispersal corridors, which has been demonstrated in many 
other species (Yanes et al., 1995). Furthermore, on both 
sides of the railway, the construction activities produced 
many sand dunes that are suitable burrowing habitat. 
These dunes might attract lizards toward the railway 
and hence increase the exchanges between the two sides 
(Figure 1). 

The temporal scale may be too small for any 
barrier effect to accumulate. It is possible that 
the railway does have a barrier effect but the 
accumulation of this effect is still below a detectable 
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level at the present time. Nevertheless, microsatellite 
DNA data and assignment tests are very sensitive, 
and are able to detect population differences 
generated in a short period of time (Pritchard et al., 
2000; Piry et al., 2004; Waser and Hadfield, 2011). For 
example, using microsatellite DNA data and assignment 
tests, Crosby et al. (2007) found significant differences 
between populations of woodfrog (Rana sylvatica) on 
both sides of Highway 401 in southern Ontario, which 
was completed 37 years prior to sampling time. In another 
study, Lesbarréres et al. (2006) found the populations 
(Rana dalmatina) separated by a highway, which has been 
open for only 21 years, had significantly more genetic 
structure (F’5;= 0.238) than the populations without a 
highway separation (F's; = 0.022) in western France. 
Considering the high sensitivity of the Plateau ecosystem, 
it is imperative to continuously monitor any changes over 
long periods of time. 

The small spatial scale may also contribute to our 
inability to detect any barrier effect of the railway. All 
of our sampling sites along the railway are within an 
approximately 20 km stretch. We, however, did detect 
differences between the railway sites and the southeastern 
sites, despite the small gap (20-25 km) between the 
two sampling areas. Other fine-scale studies of reptiles 
found significant population structure at a distance of less 
than 500 m (Sphenodon punctatus; Moore et al., 2008). 
Dileo et al. (2010) found significant population sub- 
structure of foxsnake (Mintonius gloydi) within a distance 
of 5 km. More research on the dispersal patterns of these 
lizards is needed to fully understanding the impact of the 
railway. 

Sampling bias could have masked the existing 
population differentiation. Within the railway sampling 
area in this study, many suitable habitat patches are 
nearby the railway (Figure 1), and such patches become 
fewer farther away from the railway. Consequently, many 
samples were collected fairly close to the railway from 
both sides. The short distance between sites from different 
sides of the railway may make it difficult to detect any 
barrier effect unless the effect is very strong. 


5. Conclusion 


We did not detect any significant barrier effect of 
the Qinghai-Tibet Railway on a lizard species. This 
conclusion is contradictory to the findings from studies 
on large mammal species, such as antelopes and gazelles. 
The railway structure itself may not impede lizard 
movement in a significant way, and/or the existing bridge 
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underpasses or culverts may facilitate crossing. On the 
other hand, the railway may in fact have significant barrier 
effect, but temporal scale is too small for variations to 
accumulate. Therefore, continuing the monitoring is 
important. 

The Qinghai-Tibet Plateau maintains a unique and 
fragile ecosystem, which is sensitive to disturbance. With 
the increased economic development in the region, we 
expect the traffic volume will be increased significantly, 
and more roads will likely be constructed in the near 
future. Carefully evaluating and monitoring the impact of 
such habitat alteration is crucial to protect this ecosystem. 
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Appendix I Sample site and sample size (N) information. 


Latitude (°N) Longitude (E) Group assignment Sites N Latitude (N) Longitude (E) Group assignment 


Sites N 

1 3 36.65896 95.06916 Western 25 1 36.6209 95.03953 Eastern 

2 4 36.64955 95.06145 Western 26 5 36.61631 95.03472 Eastern 

3 4 36.64612 95.0586 Western 27 1 36.60283 95.02392 Eastern 

4 3 36.64101 95.05428 Western 28 8 36.6011 95.02233 Eastern 

5 1 36.6363 95.05052 Western 29 3 36.58434 95.01044 Eastern 

6 2 36.63604 95.05025 Western 30 N 36.5827 95.0096 Eastern 

To 10 36.62437 95.03752 Western 31 7 36.56934 95.00228 Eastern 

$o all 2202009 Western 32 15 36.56903 95.0039 Eastern 

2o 2 36:00899 23:028 Western 3 2 36.5404 94.98188 Eastern 
W j cate 230121 Western 34 16 36.52213 94.96795 Eastern 
aa - : ve a A ree ifera 35 1 36.35925 95.047 Southeastern 
13 1 36.5859 95.01216 Western 36 5 36.32434 95.26437 Southeastern 
14 6 36.56981 95.00185 Western 37 3 36.33952 95.36497 Southeastern 
15 6 36.5582 94.9942 Western 38 1 36.35709 95.45422 Southeastern 
16 1 36.55534 94.99097 Western 39 1 36.37288 95.61159 Southeastern 
17 7 36.53059 94.972 Western 40 4 36.38054 95.745 13 Southeastern 
18 1 36.52745 94.9709 Western 41 1 36.37277 95.79305 Southeastern 
19 8 36.52617 94.96852 Western 42 11 36.38325 95.98775 Southeastern 
20 2 36.67026 95.07748 Eastern 43 10 36.38324 96.13561 Southeastern 
21 1 36.66687 95.08364 Eastern 44 14 36.37457 96.19838 Southeastern 
22 2 36.66637 95.0765 Eastern 45 3 36.37812 96.2617 Southeastern 
23 2 36.65367 95.06547 Eastern 46 1 36.38238 96.38237 Southeastern 
24 5 36.63562 95.05228 Eastern 47 1 36.38136 96.40668 Southeastern 


